Variational Perturbation Theory for Markov Processes 



Hagen Kleinert, Axel Pelster, and Mihai V. Putz 
Institut fur Theoretische Physik, Freie Universitat Berlin, Arnimallee 14, D-14195 Berlin, Germany 
E-mails: kleinert@physik.fu-berlin.de, pelster@physik.fu-berlin.de, putz@physik.fu-berlin.dc 

(Dated: February 1, 2008) 

We develop a convergent variational perturbation theory for conditional probability densities of 
Markov processes. The power of the theory is illustrated by applying it to the diffusion of a particle 
in an anharmonic potential. 

PACS numbers: 02.50.-r,05.10.Gg 

CN ■ 

O ' I- INTRODUCTION 

o : 

Variational perturbation theory Q transforms divergent perturbation expansions into convergent ones with the con- 
i-C . vergence extending to infinitely strong couplings. The theory has been developed for the path integral representation 
of the free energy and the density matrix in quantum statistics, and has been tested for many systems, in particular 
the anharmonic oscillator and the hydrogen atom without and with a homogeneous magnetic field Jl|, ||, ||, [|, [|, |[ . 
The procedure is based on approximating a potential by a local trial oscillator whose frequency is optimized order 
by order for each set of external end points. Recently, variational perturbation theory has also been successfully 
extended to statistical field theory to calculate highly accurate critical properties of second-order phase transitions || . 

> ' . . ,, 

qq . In this paper we develop a similarly convergent variational perturbation theory for the path integral representation of 

£ — ' the conditional probability density of Markov processes. In close analogy with the previous method we approximate a 

given stochastic process by a trial Brownian motion with a linear drift coefficient, and optimize the damping constant. 

We illustrate the procedure by calculating the time dependence of the conditional probability density for a nonlinear 

stochastic model. After some introductory remarks on Markov processes in Section II, the path integral for the 

conditional probability density is treated perturbatively in Section III, and evaluated variationally in Section IV. 

II. MARKOV PROCESSES 

We start by summarizing the basic properties of Markov processes ||, [l(], [ll], [l^] needed in the sequel. 

£h ! 

O . 

J"^ ' A. Fokker-Planck Equation 

> : 
• i— i , 

A Markov process for a single stochastic variable x is characterized by the conditional probability density P(xb tb\x a t a ) 
that the event Xb is realized at the time tt> once the event x a has taken place at time t a - It has the initial condition 

a ; 

P(x b t a \x a t a ) = 5(x b ~ X a ) (1) 

and obeys, for additive noise, the Fokker-Planck equation 

Q 

— P(xb tb\x a t a ) = L F p(xb) P{xb tb\x a t a ) , (2) 
Otb 



with the Fokker-Planck operator 



2 



L FP (x b )* = [ K ( x b)*] + D q^2* i (3) 

where K{x) and D denote the drift and diffusion coefficient, respectively. An important example is provided by the 
overdamped motion of a Brownian particle with mass M and friction constant k in an external potential V(x). In 
this case, the drift coefficient reads 



(4) 



2 



and the diffusion coefficient is proportional to the temperature T via Einstein's relation 

D = Mnk B T. 



(5) 



Since the spatial derivatives in the Fokker-Planck equation (0), (||) are all on the left-hand side, they guarantee the 
probability conservation 



d_ 

dt b 



+ 00 



dx b P{x b t b \x a t a ) = G, (6) 
such that the normalization integral, which is unity at the initial time t b = t a due to (|l]) , remains so for all times: 

dx b P{x b t b \x a t a ) = 1 ■ (7) 
In the long-time limit t b — > 00, the conditional probability density P(x b t b \x a t a ) becomes stationary 



lim P(x b t b \x a t a ) = P st (x b ), 

th — >oo 



D 



P s t(x b ) = Nexp 

where the normalization constant N follows from ((?]) and (g) 

+00 p ^ r x b 



N 



dx b exp 



D 



dx K(x) 



dx K{x) 



(8) 



where P s t(x b ) denotes the time-independent solution of the Fokker-Planck equation (0), defined by 

L FP {x b )P st (x b ) =0. (9) 
By applying the Fokker-Planck operator (^), we verify this solution is 



(10) 



(11) 



B. Path Integral 

The solution of the Fokker-Planck initial value problem has the path integral representation 

-A\x] ( 12 ) 



rx(t b ) = x b 

P{x b t b \x a t a ) = I Vxe~ 

J x(t a ) = X a 



with the generalized Onsager-Machlup functional 



1 



Ax] = T7^ dt[x(t) - K(x(t))} 2 + - dtK'(x(t)) 



*DJu- " *K ' (13) 

where all paths x(t) contribute which connect the spacetime points (x a ,t a ) and (x b ,t b ). The extra term in ([l3]) is 
needed since the path integral (|l2|), ( |l3| ) is by definition symmetrically ordered in the product of x(t) and K(x(t)), 
corresponding to a midpoint discretization Q : 

N + l 

1 



N 


+00 1 


n 

n=l 


J dx„ > 
-00 ) 




N + l r 


{- 


— y 

AD ^ 

n=l L 



AttDe 



Xn X n — \ ( X n ~t- X n —\ 

— K 



N+l 



/ / X n "T X n — 1 



Making use of the stationary solution (fL0"|), the path integral (O), (O) can be factorized as 



P(x fc t b |a; Q t Q ) = \ rJ t , Xb . P(x b t b \x a t a ) 

V ^st(Xa) 



(14) 



(15) 
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with the remaining path integral 



P{x b t b \x a t a ) 



T>x exp \ — — — I dt 



(t a )=x a 



ZD J t 



~x 2 (t) + ~K 2 (x(t))+DK'(x(t)) 



This coincides with the quantum-statistical imaginary-time evolution amplitude 



(x b hf3\x a 0) 



x{H@)=X} J 



(P)=x a 



T>x exp ■ 



up 



dr 



M 



x 2 (t) + U(x(t)) 



in which we identify 



(16) 



(17) 



t a = , t b = Hj3 , 



D = 



2M 



and set 



U(x) = ™K 2 (x) + ^K>(x). 



(18) 



(19) 



To the path integral ([17|), we can apply directly the known variational perturbation theory j]J, which will lead us in 
the present context to a solution of the Fokker-Planck initial value problem (0)-(^). 



C. Brownian Motion 



A solvable trial path integral is provided by the Brownian motion with a linear drift coefficient 

K(x) — —KX , 

where the stationary solution (|To|), ( O ) reads 



P K ,Bt(x) = 



2ttD 



exp 



(-— x 2 ) 
V 2D J 



The conditional probability density factorizes therefore according to (|l5|) as 

P K (xbtb\x a t a ) = exp - (tb - t a ) - — (x 2 - x 2 a ) P K (xbtb\x a t a ) . 



and the remaining path integral is simply 



P K (xb tb\x a t a ) 



x(t b )=x b 



x(t a )—x a 



dt 



(20) 



(21) 



(22) 



(23) 



It describes a quantum-statistical harmonic oscillator with the potential U(x) = Mk 2 x 2 /2. Inserting the imaginary- 
time evolution amplitude of the harmonic oscillator JL3| and taking into account the identification (|l8|) , we obtain 



P K (x b t b \x a t a ) 



exp 



4ttD sinh n(t b — t a ) ( 4D sinh 



— \(x 2 a + xl) cosh K(t b - t a ) - 2x a x b ] 

lK(t b - t a ) J 



(24) 



The resulting conditional probability density of the Brownian motion (|20|) follows from Q22j), (|24|) and leads to the 
well-known Gaussian distribution 



P K (x b t b \x a t a ) = , / exp ■ 

y'2TTa I (x a ,t a ;t b ) 



[x b - x(xg,t a ;t b )Y 
2a 2 (x a ,t a ;t b ) 



(25) 



with the average point 



(26) 
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and the width 



<T(x a ,t a ;t b ) = ^ [l - e -a«(*»-t.)] . (27) 

It can easily be verified that the conditional probability density (|25|)-(|27|) satisfies the initial condition (|l|), and obeys 
the Fokker-Planck equation associated with the drift coefficient pG) : 

d d d 2 

-^-P K (x b t b \x a t a ) = — {nx b P K (x b t b \x a t a )} + D^-^P K (x b t b \x a t a ) . (28) 

In general the drift coefficient K(x) of a Markov process is a nonlinear function in x, and the corresponding conditional 
probability density P(x b t b \x a t a ) cannot be calculated exactly. We must then resort to approximation methods, and 
we want to show in this paper that variational perturbation theory is a very efficient one. For this we first need an 
ordinary perturbation expansion which will be derived in the next section for a typical example. 



III. PERTURBATION THEORY 



Consider the nonlinear drift coefficient 



K(x) 



-kx — gx 



(29) 



with a coupling constant g. Such a stochastic model is useful, for instance, to describe the statistical properties of 
laser light near the threshold in semiclassical laser theory JlJ, Q. In this case, the stochastic variable x is identified 
with the electric field. The parameter k is proportional to the difference between the pump parameter a and its 
threshold value (Jthr- The coupling constant g describes the interaction between light and matter within the dipole 
approximation, and the diffusion constant D characterizes the spontaneous emission of radiation. For such a stochastic 
process, the stationary solution ([To|), (O) reads 



V K 



exp 


' 1 U 2 




9 i\ 


2D \4g ' 


f- KX 2 - 


~* ). 



1 4 {m) 



(30) 



where K v (z) denotes a modified Bessel function [lq ]. The path integral for the conditional probability density 
corresponding to (|12|), (|l3|) reads: 



P(x b t b \x a t a ) 



x(t b )=x b 



T>x exp 



' x(t a )—x a 

The decomposition of type (p5j) leads to 



I 



2 1 

dt[x(t) + nx(t) + gx 3 (t)] +- 



dt [K + 3gx 2 (t)] 



P(x b t b \x a t a )=exp -{% -t a ) - — (xl -x 2 a ) - — (xl -x*) P(x b t b \x a t a ) 



with the remaining path integral 



P{x b t b \x a t a ) 



x(tb) = Xh 



T>x exp < — 



x(ta)=X a 



2D J fa 



dt 



\i 2 {t) + X -K 2 x 2 (t) + g K x 4 (t) + ^g 2 x 6 (t) - ZgDx 2 {t) 



(31) 



(32) 



(33) 



For zero coupling constant g, we obtain the conditional probability density (|23J) of the Brownian motion (|20j). Ex- 
panding the exponential in powers of g, we find the first-order approximation 



P(x b t b \x a t a ) = P K (x b t b \x a t a )\l + g 



dt(x 2 (t)) K -— dt(x\t)) K 



2D J t 



(34) 



On the right-hand side, we have denoted the harmonic expectation value of any functional F [x] of the path x(t) by 



(F[x]) K = ^ 



P K (x b t b \x a t a ) Jx(t a )=x 



x(tb)=x b 



VxF[x] exp i — — — - / dt 



2D Jt 



\i 2 (t) + X -K 2 x 2 (i) 



(35) 



5 



The latter is evaluated with the help of the generating functional for the quantum-statistical harmonic oscillator 



x(tb) = Xb 



x(t a )=x a 



P K (x b t b \x a t a )[j] 
which has the explicit form p"3[ 

P K (xbtb\x a t a )\j] = P K (x b t b \x a t a )exp 



T>x exp 



1 

2D 



dt 



\i?{t) + - j(t)x(t) 



2D 



dhx c x(ti)j(ti) + 



4L> 2 



dh / dt a G(ii,t 2 )i(ti)j(*a) 



The quantity P^{x b t b \x a t a ) is the same as in Eq. (|24]), and x c \(ti) denotes the classical path 

x a sinh n(t b — t\) + x b sinh n(ti — t a ) 
x c \\t\) = . r 77 7^ • 

The quadratic term in (|3^) contains the Green function 

D[Q(t% — t%) sinh ft(tb — ti) sinh n(t2 — t a ) + 0(^2 — ti) sinh — t%) sinh — i a )] 



G{h,t 2 ) 



k sinh fc(tft — t a ) 



(36) 



(37) 



(38) 



(39) 



We evaluate harmonic expectation values of polynomials in x arising from the generating functional ( p7\ ) according 
to a slight generalization of the standard Wick theorem jl6| [lTj . The generalization is required by the presence of the 
linear term in (|37|). The evaluation is most economically done in a recursive procedure which we illustrate with the 
harmonic expectation value 



(x n (n)x m (r 2 )) f 



(40) 



(i) Contracting x(t\) with x n 1 (ti) and x rn (T2) leads to a Green function G(t\,ti) and G(ti,T2) with multi- 
plicity Ti—l and m, respectively. The rest of the factors remains inside the expectation symbol, leading to 
(x n - 2 {n)x m (T 2 )) K and {x n ~ 1 {Tx)x m ~ 1 (T2)) K . 

(ii) If n > 1, we extract one x(t\) from the expectation value giving x c \(ti) multiplied by (x n ~ 1 (Ti)x m (T2)) K . 
(hi) Add the terms (i) and (ii). 

(iv) Repeat the previous steps until only products of expectation values (x{t\)) k = x c \(t\) remain. 
With the help of this recursive procedure, the first-order harmonic expectation value (a; 4 (ri)) K is reduced to 

(x 4 (n)) K = x d (n) (x 3 (n)) K + 3G(r 1 ,T 1 )(x 2 (r 1 )) ti . (41) 

Similarly, we find 

(x 3 (n)) K = x c1 (t 1 )(x 2 (t 1 )) k + 2G{t 1 ,t 1 )x c1 (t 1 ), (42) 

and 

(x 2 (r 1 )) K = 2 ; c 2 1 (r 1 ) + G(r 1 ,r 1 ). (43) 
Combining Eqs. we obtain in first order 

(x\ n )) K = 4(n) +6x 2 cl (n)G(n, n) + 3 G 2 (n,Ti). (44) 

The contractions can be represented graphically by Feynman diagrams with the following rules. Vertices represent 
the integrations over t 



dt. 



dt. 



(45) 



a line denotes the Green function 



1 2 = G(il,i 2 ), 



(46) 
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and a line ending with a cross represents the classical path 

x 1 = x c i(ti). 



(47) 



Inserting the harmonic expectation values ( f43| ) and (|44|) into the perturbation expansion (|34|) leads to the first-order 
diagrams 



P(x b t b \x a t a ) = P K {x b t b \x a t a ) I 1+g 



K 

2D 



■(•48) 



Evaluating the diagrams, the conditional probability density for the drift coefficient ([29|) following from (32) and ( |4S| ) 
becomes to the first order in the coupling constant g 



(x b t b \x a t a ) = P K {x b t b \x a t a ) jl + g c$(t) + c^q{t)xI + (r)x a x b + (r)a 



+ c< 40 ( T ) X a + C 41 ( T ) x a x b + c 42 ( T ) X a x b + c 43 ( T ) x a x b + c 44 ( T ) X b 



„(1 



(!)/ 



where the expansion coefficients c|^(t) are the following functions of the dimensionless variable r = n(t b — t a ) (see 



(49) 



Fig. @): 



3D [I + 4e 


- 2T (I-2r)-e- 4T (5 + 4r)] 






4K 2 (I-e- 2T ) 2 




3 [e- 2r (4r 


-5) + 4e- 4r (l + 2r) + e- 6r ] 


c£>(r 




2n{l-e- 2T f 


3 [e" r (2- 


t) + 2e- 3r (I - At) - e - 5T (3r + 


4)] 




k(1 -e^) 3 


5 


2e" 2r + 3e 


- 4r (l-4r)-6e- 6T + e- 8T 






4D(1 - e^) 4 




-e~ T + 3e 


" 3t (4t - 3) + 3e" 5r (3 + 4r) + e 


-7r 




2D(1 - e- 2 ^) 4 




3 [e- 2r (3 - 


- 2r) - 8re- 4r - e - 6T (3 + 2r)] 






2 J D(I-e- 2 ^) 4 




-I + 6e~ 2r - 3e- 4r (I + 4r) - 2e~ 6r 






4D(1 - e" 2 ^) 4 





c (1 Vr 



c (1) fr 



c 44 l T 



It is easy to verify that the conditional probability density (f49|)-(|50|) 
(i) obeys the corresponding Fokker-Planck equation following from (|^) 

d d d 2 

— P(x b t b \x a t a ) = ^— [{KXb + gX b ) P{x b t b \x a t a )] + D-^-^P(x b t b \x a t a ) 

with the initial condition (liT), because of 



(50) 



and (M 



(51) 



(0) = cg } (0) = <&> (0) = c&> (0) = c£ (0) = 4 1 ' (0) = (0) = , c# (0) = - C « (0) = ^ . (52) 

(ii) is normalized according to (R) for all times ib. 

(iii) approaches the stationary solution |3(]) in the long-time limit (^) . 

All these properties of P{x b t b \x a t a ) are satisfied to first order in the coupling constant g. As such the perturbative 
result (p9|)~(p0|) is only applicable for small values of the coupling constant g. This limitation is removed by a 
variational evaluation which enables us to find an approximate conditional probability density for all values of the 
coupling constant g. 
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4i>w 

0.15 



-0.15 




8 10 



4o>M 




1 2 3 4 5 6 




C «M 2 4 6 8 10 




1 2 3 4 5 6 





FIG. 1: Temporal behavior of the expansion coefficients ( |5o| ) of the conditional probability density ( jis)) ) as a function of 
r = tt(tb — t a ) for D = k = 1. 



IV. VARIATIONAL PERTURBATION THEORY 



Variational perturbation theory in quantum statistics approximates a given potential by adding and subtracting a 
local trial harmonic oscillator whose frequency is then systematically optimized . Here we transfer this procedure to 
Markov processes and modify the nonlinear drift coefficient (^9|) by adding and subtracting a trial Brownian motion 
with a yet unknown damping constant K: 

K{x) = -Kx - g {^—^ x + x3 J ■ ( 53 ) 

After this, we treat the combined second term as being small of the order of the coupling constant g. The result 
is obtained most simply by replacing the damping constant n in the original drift coefficient ( |29| ) according to the 
substitution rule (compare with Chap. 5 in Ref. [Q) 



where we have introduced the abbreviation 



K(l + gr) , 



k-K 



gK 



Writing the conditional probability density (|4^) as 

P{x b t b \x a t a ) = exp [W{x b t b ;x a t a )] 



(54) 



(55) 



(56) 



8 



we apply the substitution rule (54) to the cumulant expansion W(xb H; x a t a ), and reexpand up to the first order in 
the coupling constant g. Afterwards the abbreviation r is reexpressed in terms of n and K via (|55|). This leads to 

(x b t b - x a t a -K) = {c$ (r) + c$ (r)x 2 a + 4? (r)^ + 4? (r)x 2 b + g [c$ (r) + c&> (r)^ + <&> ( T ) Xa2;b 

+4 2 V)^+4oV)-^+4iV)^ 

where the zeroth-order expansion coefficients 



4^ (r)a£ X 2 + c$ (r)x a zf + c# (t)^] } , (57) 



4°o } (r) 



1, ^_ 

2 n 27rD(l- 



Ke 2t (k — if )re 

2L>(1 - e- 2r ) + L>(1 - e" 



1 

-2t 



-2t 



„(0) 



„(0) 
!., 22 



(r) = 



(K- J ftT)T(l + e- 2T )e- 



D(l 



-2r\ 



(k — K)re 



-2t\2 



-2t 



2D{l-e- 2T ) D(l-e- 2 ^) 2 



(58) 



and the first-order expansion coefficients ( |50| ) are functions of the dimensionless variable r = K(t\, — t a ). Note that 
using (|50|) in (b7|) necessitates to substitute k by if. 



We now remove the dependence of the cumulant expansion (57) from the artificially introduced trial damping constant 
K. According to the principle of minimal sensitivity [ |j"8| , we minimize its influence on W^- 1 ' (xi, tb\x a t a ; K) by searching 
for local extrema of W^ 1 ' (x b tb\x a t a ; K) with respect to K, i.e. from the condition 



dW^(x b tb;x a t a -K) 



dK 



(59) 



K=K^(x b t b ;x a t a ) 



It may happen that this equation is not solvable within a certain region of the parameters XbAb, x a ,t a . In this case 
we look for turning points instead in accordance with the principle of minimal sensitivity [[l], [l?]], i.e. we determine 
the variational parameter K^{xbtb\x a t a ) from solving 



d 2 W^{x b t b ;x a t a ;K) 



dK 2 



= 0. 



(60) 



K=K( 1 )(x b t b -x a t a ) 

The solution K^\xb h; x a t a ) from (59) or ( |60| ) yields the variational result 

exp [W (x b t b \X a t^K^ixbt^Xa t a ))] 



P(x b t b \x a t a ) 



dxb exp 



IT' 



{^X}) £fr, X a £ a , ^ {xb ^bi ^a) 



(61) 



for the conditional probability density. Note that variational perturbation theory does not preserve the normalization 
of the conditional probability density. Although the perturbative result (fj^-([50|) is still normalized in the usual 
sense (Q) to first order in the coupling constant g, this normalization is spoilt by choosing an x^-dependent damping 
constant K^\xb W, x a t a ). Thus we have to normalize the variational conditional probability density according to 
( |6l| ) at the end (compare the similar situation for the variational ground-state wave function in Refs. Jl?], |l9[]). 



We have applied variational perturbation theory to analyze the nonlinear stochastic model (|2S|) with D = 1, x a = t a = 
below and above the laser threshold, i.e. k = 1 and k = — 1, for small and strong coupling g — 0.1 and g — 10. The 
results for the conditional probability density P(xb t a \x a t a ) are plotted in Figs. and |^. On the scale of the figures, 
they show no significant deviation from numerical solutions of the corresponding Fokker-Planck equation, some minor 
deviations only occur above the laser threshold k = — 1 for g = 0.1 p3[ - Both figures illustrate how the densities, 
originally peaked at the origin, turn into their stationary solutions in the long-time limit. The stationary solution 
reveals below the threshold one extremum in Fig. whereas above the laser threshold in Fig. [| two extrema occur. 
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V. CONCLUSION 

We have presented the lowest-order variational calculation for the conditional probability density of the nonlinear 
drift coefficient (p9|). By going to higher orders in variational perturbation theory, it is straightforward to increase 
systematically the accuracy to any desired degree |IJ . 
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Note Added in Proof 

Meantime we became aware of the alternative approach j2^] to variationally analyze the nonlinear model (^9j). Above 
the laser threshold k = — 1, that method yields for all times if, a unique solution of the extremal condition corresponding 
to (^9|). However, the resulting conditional probability density shows for larger times if, significant deviations from 
our and from numerical solutions of the Fokker-Planck equation. 

P(x b t b \x a t a ) P(x b t b \x a t a ) 




FIG. 2: Conditional probability density P(x b t b \x a t a ) below the laser threshold at k = 1, D = 1, and x a — t a — 0. For the 
coupling constants g = 0.1 in a) and g = 10 in b) the distribution is shown for the times t b = 0.05, 0.1, 0.2, 0.5, 1 from the top 
to the bottom at the origin, respectively. 

P(x b t b \x a t a ) P(x b t b \x a t a ) 




FIG. 3: Conditional probability density P(x b t b \x a t a ) above the laser threshold at k = — 1, D = 1, and x a = t a = 0. For the 
coupling constants g = 0.1 in a) and g — 10 in b) the distribution is shown for the times t b = 0.1,0.2,0.5,1,1.5,2,3,4 and 
t b = 0.05, 0.1, 0.2, 0.5 from the top to the bottom at the origin, respectively. 
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